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ABSTRACT 

We describe a technique that finds orbits through the Galaxy that are consistent with 
measurements of a tidal stream, taking into account the extent that tidal streams 
do not precisely delineate orbits. We show that if accurate line-of-sight velocities are 
measured along a well defined stream, the technique recovers the underlying orbit 
through the Galaxy and predicts the distances and proper motions along the stream 
to high precision. As the error bars on the location and velocities of the stream grow, 
the technique is able to find more and more orbits that are consistent with the data 
and the uncertainties in the predicted distances and proper motions increase. With 
radial- velocity data along a stream ~ 40° long and <0.3° wide on the sky accurate to 
~ lkms - the precisions of the distances and tangential velocities along the stream 
are 4 percent and 5kms _1 , respectively. The technique can be used to diagnose the 
Galactic potential: if circular-speed curve is actually flat, both a Keplerian potential 
and <3>(r) oc r are readily excluded. Given the correct radial density profile for the dark 
halo, the halo's mass can be determined to a precision of 5 percent. 
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1 INTRODUCTION 



Deep optical surveys of the Milky Way and other 
Local- Gr oup galaxies have uncovered numerous stellar 
streams (lOdenkirchen et al.l 2001 : FMaiewsk i et al ] 120041 : 



iBelokurov et all 120061 : llbata et all I2007T ) The Leiden- 



Argentin e- Bonn survey of th e Galaxy in the 21 cm line of hy- 
drogen (|Kalberla et al.ll2005h contains many similar streams. 
In all probability both stellar and gaseous streams have been 
tidally torn from orbiting bodies, and as suc h delineate the 
orbit s of those bodies around the Galaxy (I Johnston et aD 
Il996l : lOdenkirchen et al.ll20Q3l : IChoi et al.ll200/t ). Newton's 
laws of motion severely constrain the readily observable 
quantities along an orbit in the sky, namely the sequence of 
positions on the sky [l(u),b(u)] and the corresponding line- 
of-sight velocities, v\\(u), where u is a parameter that varies 
monotonicall y along the stream (|jin Sz Lvnden-Belll 120071 : 
Binn evl 120081 . hereafter Paper I). In fact, if the observables 
are known to reasonable accuracy, data for a single stream 
can strongly constrain the Galaxy's gravitational potential, 
and once the potential is known, the distance and proper 
motion at each point on the stream can be predicted with 
an accuracy that far exceeds anything likely to be possible 
by conventional astrometry (Paper I). 

From the work of Paper I it emerges that the major lim- 
itation on the diagnostic power of streams is that streams d o 
not precisely delineate individual orbits (|Choi et al.l l20Q7h . 
This paper is devoted to exploring the extent to which this 



limitation can be overcome. In Section [2] we illustrate the 
extent of the problem, in Section [3] we introduce significant 
improvements to the methodology of Paper I and use these 
to identify orbits that are consistent with a given body of 
data. In Section [4] we test this approach. In Section [5] we 
examine our ability to correctly diagnose the Galactic po- 
tential. Section [6] sums up and discusses directions for future 
work. 

Except where stated otherwise, orbits and reconstruc- 
tions are calculated using the G alactic potential of Model 
II in iBinnev &; Tremaind (|2008l ) , which is a slightly mod- 
ified version of a halo-do minated potential described by 
iDehnen fe Binnevl ([l998a). We take the distance to the 
Galac tic centre to be 8kpc and from Reid Brunth aler 
(|2005l ) (for V and W) and IDehnen fc Binnevl (|l998bh we 
take the velocity of the Sun in the Galactic rest frame to be 
(U,V, W) = (10.0,241.0,7.6) kms" 1 . 



2 THE PROBLEM 

The full curves in Fig.[T]show an orbit superficially similar to 
that underlying the Orphan Stream (Bel okurov et aLl feoOT) 
from two viewing locations - the position of the Sun and a 
position 120° further round the solar circle. Also shown in 
each projection are the locations of particles tidally stripped 
from a self- gravitating N-body model of a cluster launched 
onto the given orbit. Clearly the particles provide a very 
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Figure 1. Full (red) lines: the orbit of a progenitor of an Orphan-like stream. Broken (green/blue) lines: orbits of a star now seen at either 
end of the tidal tail. Points: particles tidally stripped from an N-body model of the Orphan-like progenitor. Upper panels: distributions on 
the sky; lower panels: line-of-sight velocities. The N-body model had 60 000 particles set up as a King model with W = 2, ro = 13.66 pc 
and Mo = 9381MQ, on the orbit detailed in Tabled and evolved for 9.43 Gyr. The particles were advanced in time by the "FVFPS" 
tree code of Londrillo et al. (2003). 



useful guide to the orbit of the cluster, but they do not pre- 
cisely delineate it. Moreover, the relationship of the orbit 
to the stream depends on viewing angle. The line-of-sight 
velocities of stream particles have a similar relationship to 
the orbit's line-of-sight velocity. Hence even with perfectly 
error-free data, the orbit we seek will not coincide with the 
stream, and before we can fully exploit the dynamical po- 
tential of streams, we have to understand how to infer the 
location of an underlying orbit from measurements of the 
stream. 

With what precision could the track of an orbit be spec- 
ified from the positions of the particles in Fig. [TJ First we 
need to be clear that any orbit will do. Generally it will be 
convenient to use the orbit that passes through some point 
that lies near the centre of the observed stream, both on 
the sky and in line-of-sight velocity. In some circumstances 
this orbit will closely approximate the orbit of the centre of 
mass of the stream's progenitor, but there is no requirement 
that this is so. Since a central point on the orbit is chosen 
at will, this point is associated with vanishing error bars. 
As one moves away from this point, either up or down the 
stream, it becomes more uncertain where the chosen orbit 
lies, and the size of our error bars must increase. Hence the 
region of (/, b, v\\) space to which the orbit is confined by the 



observations is widest at its extremities and shrinks, usually 
to a point, at its centre. We call this the "bow-tie region". 

In the top right panel of Fig. [T]the leading and trailing 
streams are not offset for most of the span so we may assume 
that the orbit through the stream's centre runs right down 
the centre of the stream. In the lower panels the stream has 
a kink at the progenitor and our best guess is that the orbit 
through the point where the two halves of the tail touch runs 
near the lower edge of the left-hand half and near the upper 
edge of the right-hand half. In every case the error bars on 
the location of the orbit grow from zero at the centre to 
roughly half width of the stream at its ends. Quantitatively, 
the largest error is then 0.15° for the top- left panel, ~ 0.25° 
for the top-right panel, and ~ 5 km s _1 in the bottom panels. 



3 IDENTIFYING DYNAMICAL ORBITS 

Paper I showed that given an orbit's projection onto the 
sky [l(u),b(u)] and the corresponding line-of-sight velocities 
v\\(u), the remaining phase space coordinates can be recov- 
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ered by solving the differential equations 




ds _ dt 
du " du 

Here u is distance on the sky down the projected orbit, s(u) 
is three-dimensional distance to the orbit, t(u) is the time 
at which the orbiting body reaches the given point on the 
orbit and F\\ is the component of the Galaxy's gravitational 
field along the line of sight. The reference frame used is 
the inert ial frame in which the Galactic centre is at rest; 
consequently the velocities v\\ are obtained by subtracting 
the projection of the Sun's motion along the given line of 
sight from the measured heliocentric velocities. 

If the input data used to solve equations (pQ) are not de- 
rived from an orbit in the same force field as is used to derive 
F|| , the reconstructed phase-space coordinates will not sat- 
isfy the equations of motion. Paper I observed that violation 
of the equations of motions might cause the reconstructed 
solution to violate energy conservation, and therefore used 
rms energy variation down the track as a diagnostic for the 
quality of a solution. However, energy conservation is neces- 
sary but not sufficient to qualify a track as being an orbit. 
Here we construct a diagnostic quantity from residual er- 
rors in the equations of motion themselves, since orbits are 
defined to be solutions of these equations. 

We first derive the equations of motion. In the coordi- 
nates (s,6, /), the canonically conjugate momenta are 



Ps = s 
Pb 



s 6, 



2 2, ; 
pi — s cos bl. 

The Hamiltonian is therefore 
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and the equations of motion are 
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As in Paper I, when solving equations (pQ) extensive use 
is made of cubic-spline fits to the data. In the examples 
presented in Paper I natural splines were used in order to 
avoid specifying the gradient of the data at its end points. 
Significantly improved numerical accuracy can be achieved 
by taking the trouble to specify these gradients explicitly. 
Given input data, we estimate the quantity dl/db at the 
end points by fitting a quadratic curve through the first 
three and last three points, dl/db is then computed at the 
location of the middle point of each set, and the very first 
and very last points are considered 'used' and thrown away. 
This quantity is then used in the geometric relation 



du 
db 



■ ±W 1 + cos 2 b 



(5) 



to compute db/du at the ends of the track. The sign am- 
biguity is resolved by inspection of the directionality of the 
input data. We then use the geometric relation 



dl_ 

du 



zb sec I 



1 



(6) 



to obtain dl/du at the ends of the track, where the sign 
ambiguity is resolved in the same way. We are now able to 
fit cubic splines through the input tracks, with the slopes at 
the end of the l(u) and b(u) tracks given as above, but at 
this stage the track of v\\ (u) is fitted with a natural spline. 
The reconstruction equations (pQ) are now solved for t(u), 
which is then fitted by a cubic spline, with the slopes at the 
ends given explicitly by JTJ). 

We can now compute l(t) and b(t) and fit splines to 
them, with the slopes at the ends computed from dl/du and 
db/du by the chain rule. The momenta J2} are now calcu- 
lated explicitly, using the derivatives of the splines l(t),b(t) 
in place of /, b. The slopes at the endpoints, dv\\/du, can 
now be calculated from @ and dt/du; the v\\(u) spline is 
refitted using these boundary conditions, the reconstruction 
repeated, and the momenta recalculated. 

The left- and right-hand sides of the equations of motion 
J4j) are calculated explicitly. For each equation of motion we 
define a residual 



R(t) =P\hs(t) -Prhs(t). 



(7) 



These residuals are used to compute, for each equation of 
motion, the diagnostic quantity 



D = lo gl 



His 



(8) 



where the residuals have been normalised by the mean- 
square acceleration and the times t\ and £2 correspond to 
the fifth and fifth- from- last input data points; the residual 
errors from the end regions, < t < t\ and £2 < t < £ ma x, 
tend to dominate the integrated quantity and are not easily 
reduced by modifying the input; they are therefore excluded. 
The largest of the three values for D is used as the diagnostic 
quantity for that particular input. 



3.1 Parametrising tracks 

Our strategy for identifying a stream's underlying orbit is 
to compute the diagnostic D (eq. [8} for a large number of 
candidate tracks, and to find which candidates yield values 
of D consistent with their being dynamical orbits. 

We start by specifying a baseline track across the sky 
[lh(u), bh(u')]j where u is a parameter that increases mono- 
tonically down the track from —1 to 1. Similarly, we specify 
associated baseline line-of-sight velocities v\\\>(u'). The base- 
line track is required to pass through the error bars of every 
data point. 

All candidate tracks should be smooth because orbits 
are. We satisfy this condition by expressing the difference be- 
tween the baseline track and a candidate track as a low-order 
polynomial in u . For streams that cover a wide range of lon- 
gitudes, the parametrisation of candidate tracks is achieved 
by slightly changing the values of b and v\\ associated with 
a given value of / from the values specified by the baseline 
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functions. That is we write 



b(u') =b h (u') + Y,bnT„(u'), 

n=0 
N 

V\\ (u) = V\\ h (u) + ^ 0>nT n (u), 



(9) 



where T n is the n th -order Chebyshev polynomial of the first 
kind and a n and b n are free parameters. These 2N parame- 
ters are coordinates for the space of tracks that we have to 
search for orbits. When a stream does not stray far from the 
Galactic plane, candidate tracks are best parametrised by 
adjusting the baseline values of b and v\\ at given longitude. 
In all examples in this paper, the series in equations J9} are 
truncated after N = 10. A larger number of terms allows 
the correction function to produce tracks that represent or- 
bits better, but makes the search procedure computationally 
more expensive. The number of terms used is a compromise 
between these considerations. 

The space of tracks is defined by the a n and b n and 
one extra parameter, the distance to the stream, so, at the 
starting point u = for the integration of equations (pQ). 

We shall henceforth denote a point in the (2N + 1)- 
dimensional space of parameters by %. Each x is associated 
with a complete specification of all six phase-space coordi- 
nates for every point on the candidate orbit: /, b and v\\ follow 
from the parametrisation and the remaining coordinates are 
obtained by solution of the differential equations (pQ). Con- 
sequently, each x corresponds to a value of D (eq. [8} that 
quantifies the extent to which the phase-space coordinates 
deviate from a dynamical orbit in the given potential. 



3.2 Searching parameter space 

Dynamical orbits are found by minimising the sum 
D'(x) = D( X )+p(x) 

where p(x) is the sum of the penalty functions: 



where 

Pi,pos - 

with 







if Ai. 



> 1 
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Here 5bi is the width in b of the bow-tie region at U. Similarly 

„ n\ / Ai ve \ if Ai vel > 1 f-\A\ 

L otherwise 

with 

_ \v\\(h) - v\\ h (h)\ 

^i,vel 



5v\\(k 



(15) 



Prior information about the distance to the stream is used 
by specifying the penalty function p s to be 



Ps 



'(3 \s - s oh \ /5s if \s - s h\ > 5s 







otherwise, 



(16) 



where 5 s is the half- width of the allowed range in the dis- 
tance so to the starting point of the integrations and sob 
is the baseline value of so- These definitions are such that 
p(x) — so long as the track lies within the region that is ex- 
pected to contain the orbit, and rises to unity, or in the case 
of p s to /3, on the boundary and then increases continuously 
as the orbit leaves the expected region. 

In practical cases the prior uncertainty in distance is 
large, and the obvious way to search for orbits is to set 5s to 
the large value that reflects this uncertainty and then set the 
algorithm described below to work. It will find candidate or- 
bits for certain distances. However, we shall see below that it 
is more instructive to search the range of possible distances 
by setting 5s to a small value such as 0.5 kpc and searching 
for orbits at each of a grid of values of sob- In this way we 
not only find possible orbits, but we show that no accept- 
able orbits exist outside a certain range of distances. In this 
procedure the logic underlying 5s is very different from that 
underlying 5b and 5v\\ . 

Since p(x) is added to the logarithm of the rms errors 
in the equations of motion and increases by of order unity at 
the edge of the bow-tie region, the algorithm effectively con- 
fines its search to the bow-tie region, where p = 0. Thus at 
this stage we do not discriminate against orbits that graze 
the edge of the bow-tie region in favour of ones that run 
along its centre. Our focus at this stage is on determin- 
ing for which distances dynamical orbits can be constructed 
that are compatible with the data. Once this has been estab- 
lished, distances that lie outside some range can be excluded 
from further consideration. 

The space of candidate tracks x is 21-dimensional, so 
an exhaustive search for minima of D' (eq. I lQf) is impracti- 
cal. Furthermore, the landscape specified by D' is complex. 
Some of this complexity is physical; the space should contain 
continua of related orbits, and ideally D' — > — oo at orbits. 
Hence deep trenches should criss-cross the space. Superim- 
posed on this physical complexity is a level of numerical 
noise arising from numerical limitations in the computation 
of D'(x)- The limitations include the use of finite step sizes 
in the solution of equations (pQ) and the subsequent evalua- 
tion of D'(x)i as we ll as the difficulty in representing a true 
orbital track with a collection of sparse input points inter- 
polated with splines. In practice numerical noise sets a lower 
limit on the returned values of D'(x)- 

On account of the complexity of the landscape that D' 
defines, "greedy" optimisation methods, which typically fol- 
low the path of steepest descent, are not effective in locating 
minima. The task effectively becomes one of global minimi- 
sation, which is a well studied problem in optimisation. 

We have used the variant of the Metropo l is "s imulated 
annealing" algorithm described in Pre ss et al . (2002), which 
uses a modified form of the downhill simplex algorithm. In 
the standard simplex algorithm, the mean of the values of 
the objective function over the vertices decreases every time 
the simplex deforms. In the Press et al. algorithm the sim- 
plex has a non- vanishing probability of deforming to a con- 
figuration in which this mean is higher than before. Conse- 
quently, the simplex has a chance of crawling uphill out of a 
local minimum. The probability that the simplex crawls up- 
hill is controlled by a "temperature" variable T: when T is 
large, uphill moves are likely, and they become vanishingly 



© 2008 RAS, MNRAS 000. [TIP 



Locating the orbits delineated by tidal streams 5 



rare as T — > 0. During annealing the value of T is gradually 
lowered from an initially high value towards zero. 

One vertex of the initial simplex is some point Xguess and 
the remaining 2N + 1 vertices are obtained by incrementing 
each coordinate of Xguess in turn by a small amount. For the 
coefficients of To this increment is approximately the size of 
the allowed half widths, 5s, 5v and 5b. Increments for coordi- 
nates representing coefficients of higher-order T n are scaled 
as 1/n. The overall size of these increments is therefore set 
by the size of the region within which we believe the global 
minimum to lie. It is important to note that in each genera- 
tion of a simplex, the increments should independently have 
equal chance of being added to or subtracted from the values 
of Xguess so that no part of the parameter space is unfairly 
undersampled. The algorithm makes tens of thousands of 
deformations of the simplex while T is linearly reduced to 
zero. 

This entire process is repeated some tens of times, after 
which we have a sample of local minima that are all obtained 
from Xguess- 

We now update Xguess to the location of the lowest of 
the minima just found and initiate a new search. The entire 
process is repeated until the value of the diagnostic func- 
tion D'(x) hrt s a floor. When this floor lies higher than the 
numerical-noise floor, the attempt to find an orbit that is 
consistent with the assumed inputs has been a failure and 
we infer that no such orbit exists. When the floor coincides 
with the numerical-noise floor, we conclude that the corre- 
sponding x specifies an orbit that is compatible with the 
inputs. An approximate value for the numerical noise floor 
for a given problem may be obtained as follows: given input 
that perfectly delineates an orbit in the potential in use, 
the value of D' returned at the correct distance is approx- 
imately the numerical-noise floor. Conclusive proof that a 
candidate track with a particular value of D' is an orbit can 
be obtained by integrating the equations of motion from the 
position and velocity of any point on the track and ensuring 
that the time integration essentially recovers the track. 

On account of the stochastic nature of the algorithm, 
an attempt to find a solution at a particular distance occa- 
sionally sticks at a higher value of D' than the underlying 
problem allows. This condition is identified by scatter in 
the values of D' reached on successive attempts and by in- 
consistency of these values with the values of D' achieved 
for nearby distances — we see from Fig. [5] that the function 
underlying the minima is smooth. When the magnitude of 
this scatter is significant, one can only confidently declare 
an attempt to find an orbit a failure if the D' achieved is 
consistently higher than the noise floor by more than the 
scatter; since the diagnostic measure D' quantifies the ex- 
tent to which a candidate track satisfies the equations of 
motion, by definition, tracks with higher D' than the noise 
floor plus scatter cannot represent orbits. 

When the observational constraints are weak, we expect 
several orbits to be compatible with them. In particular, we 
will be able to find acceptable orbits for a range of initial 
distances so. It is therefore important, for any given input, 
to run the algorithm starting from many different values of 
sob with 5s set to prevent the algorithm straying far from 
the specified sob- In this way, the full range of allowable 
distances can be mapped out, and dynamical orbits found for 
each distance in that range. In the case of significant scatter 
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Figure 2. Lines: Input track 6(7), as used for data sets PD2-PD6, 
along with upper and lower bounds as set by the penalty function 
Ppos (eq 1120 . Dots: the projection of the N-body simulation onto 
the sky, from which the input was derived. 




-150 1 1 1 1 1 1 1 1 1 1 1 1 

190 195 200 205 210 215 220 225 230 235 240 245 

1 / degrees 

Figure 3. Lines: Input track v\\(l), as used for data sets PD2- 
PD4, along with examples of upper and lower bounds as set for 
PD2 (narrow) and PD4 (wide). The bounds are enforced by the 
penalty function p ve \ (ea 1 14 J) . Dots: the projection of the N-body 
simulation onto the sky, from which the input was derived. 

about the noise-floor, the range of distances at which valid 
orbits are found is the range within which solutions yield 
values of D' smaller than the noise- floor plus scatter. 

Similar degeneracies in the parameters controlling the 
astrometry and line-of-sight velocities are less of a concern 
because if we have orbits that differ in these observables, 
we simply concentrate on the orbit that lies closest to the 
baseline track. 



4 TESTING THE METHOD 

To test this method, we used the N-body approximation to 
the Orphan Stream described in Fig. [T] as our raw data. Sets 
of points of [l,b] and [l,v\\] were selected by eye to lie down 
the middle of the stream. These sets were each fitted with a 
low-order polynomial to ensure smoothness, and these poly- 
nomials were sampled at 30 points to produce the baseline 
input data b(l) and v\\(l). To each data set we attached er- 
ror estimates 5b and 5v\\ , which through the penalty func- 
tions ppos and pvei feas 1121 and I14[) constrain the tracks that 
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Figure 4. Lines: Input track vu(l), as used for data sets PD5 
(dotted green) and PD6 (dashed blue). The unmodified input of 
PD2 (full red) is shown for comparison. Dots: the projection of 
the N-body simulation onto the sky, from which the input was 
derived. 
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Figure 5. Values of the diagnostic function D' for candidate 
orbits reconstructed from the pseudo-data set PD1. Each group 
of crosses is associated with one of 15 ranges within which the 
starting distance so was constrained to lie. For each such range 
the candidate orbit was reconstructed from 280 trial tracks, and 
each track yields a cross. In this figure the crosses are largely 
superimposed because the errors in the data set are small and 
there is little scope for tweaking the track. 



the Metropolis algorithm can try. Details of the resulting 
pseudo-data sets are given below, and are summarised in 
Table [2] 

In one case, PD1, the above baseline input data were 
replaced by those of a perfect orbit and Sb and 5v\\ set very 
narrow (6arcsec and 2 x 10 -3 kms -1 ) in order to validate 
the reconstruction algorithm. 

The uncertainty 5b(l) takes the same value for all the 
remaining pseudo-data sets because we assume that the as- 
trometry is sufficiently precise for the error in position to be 
dominated by the offset of the stream from an orbit. In all 
cases, 5b has a maximum value of 0.15° at the ends of the 
stream, falling linearly to zero at the position of the pro- 
genitor, consistent with the orbit of the progenitor seen in 
Fig. [T] Fig. [2] shows this input alongside the N-body data 
from which it was derived. 

For the pseudo-data sets PD2 and PD3, Sv\\ is set to 
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Figure 6. The upper panel shows heliocentric distance, s, versus 
galactic longitude, for the best reconstruction from Fig. [5] and 
for the true orbit from which the input was generated. The two 
curves are close to overlying. The lower panel shows tangential 
velocity, v±, for the best reconstruction and for the true orbit: 
the departure of the reconstruction from the orbit near the end- 
points is symptomatic of the problems near the endpoints that 
necessitate their excision from the diagnostic. 



a maximum at the ends of the stream, and falls linearly 
to a minimum of zero and 2km/s respectively, at the posi- 
tion of the progenitor. These examples represent the case 
in which the uncertainty in radial velocity is dominated by 
the orbital offset, and the case in which there is a significant 
contribution of random error at a level that is easily obtain- 
able with a spectrograph. In one case, PD4, 5v\\ is held fixed 
at lOkms -1 across the range. This case is representative of 
the imprecise radial- velocity information currently available 
from the SDSS for stars in distant streams. Fig.[3]shows the 
input for these data sets. 

For the pseudo-data sets PD5 and PD6, we added to the 
baseline data systematic offsets in vu (I) to mimic systematic 
errors in radial velocity. 5v\\ varies between a maximum and 
a minimum as in PD2 and PD3, with the values set to en- 
compass the (assumed known) systematic bias. Fig.|4]shows 
the input for these data sets. 

The pseudo-data set PD7 is identical to that of PD2, 
except that the number of raw [l,v\\] points was reduced 
to just three: one at either end of the N-body stream, and 
one at the location of the progenitor. A quadratic curve was 
perfectly fitted through these three points, and sampled at 
30 locations to produce the baseline v\\ (I) input. 5v\\ is set to 
the same maximum value as PD2 at the outermost points; 
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Figure 7. Upper panel: projections onto the sky of two candidate 
orbits and the N-body data from which the input was derived. 
The full curves show the candidate orbit at 46kpc from PD2, 
and the dashed curves show the candidate orbit at 43kpc from 
PD6,. Lower panel: radial velocity down the track for the same 
candidate orbits. The penalty function (eq. Ill)) forces the sky 
projection and radial velocity curves of candidate orbits to be 
consistent with the data; the equivalent plots for all other tracks 
are very similar to these examples. 



Table 1. Parameters of highlighted orbits from this paper. The 
coordinate system used is right-handed with x pointing away from 
the Galactic centre and y opposite the sense of Galactic rotation. 



Orbit position (x,y,z) /kpc velocity (x,y,z) /kms 1 

N-body Orphan (28.1, -10.0, 34.0) (-0.431, -0.179, -0.368) 
PD1 Test Orbit (35.5, 7.80, 37.8) (-0.0385, 0.272, 0.231) 



5v\\ is set to zero at the centre point. Only these three points 
are allowed to contribute to the penalty function (fl4|) in this 
pseudo-data set. 

In all of the examples, the penalty function (eq. Ill)) 
acts to constrain the candidate tracks to be consistent with 
the data. The contribution of the penalty function to D' is 
therefore zero in all examples. All candidate tracks are guar- 
anteed to be consistent with the data, even if they do not 
represent dynamical orbits. Fig. [71 provides example plots 
of (Ifb) and (l,v r ) for candidate tracks from PD2 and PD6 
along with the raw N-body data from which the baseline 
data are derived. Equivalent plots for all candidate tracks 
are similar to these. 

For PD1, PD2 and PD3, and each of 15 values of the 



Table 2. Configurations for tests of the method. 



Set 


<^||, max/km S 1 


<^||,min/kms 1 


offset i; || /kms 1 


PD1 


2 x KT 3 


2 x 10-3 





PD2 


4 








PD3 


6 


2 





PD4 


10 


10 





PD5 


6 


2 


2 


PD6 


15 


10.5 


10 


PD7 


4 









baseline distance sob, 280 optimisation attempts were made, 
each involving Metropolis annealing for 24,000 simplex de- 
formations, from an initial temperature of 0.5 dex. The start- 
ing distances were constrained by the penalty function p s 
(eq. fT6]) with f3 = 10 6 and 5s = 0.5 kpc, so the Metropo- 
lis algorithm could only explore a narrow band in so- Af- 
ter 40 optimisation attempts from a given starting point 
Xguess, the starting point was updated to the end point of 
the most successful of these optimisations, and annealing 
recommenced from a high temperature. In total 6 of these 
updates were performed. With these parameters, a search at 
a single distance completes in 12 CPU-hours on 3GHz Xeon- 
class Intel hardware. PD4, PD5 and PD6 follow almost the 
same schema, except that 48,000 simplex deformations were 
formed for each of 60 attempts at the same x> which was 
updated 6 times. A single distance in the latter case took 36 
CPU- hours to search: the computational load scales linearly 
with the number of deformations considered. 

Fig. [5] shows the results obtained with PD1, which has 
very small error bars. The diagnostic function D' has a 
smooth minimum that pinpoints the distance so — 47.4 kpc 
to the starting point with an uncertainty of ~ 0.2 kpc. 
The Metropolis optimisation can significantly reduce D' by 
tweaking the input track only when so is close to the truth. 
The depth of the minimum indicates the numerical noise 
floor for this particular problem, and no significant scatter 
is seen between successive runs. We do expect the noise floor 
to be slightly different for PD2-PD6 because both the un- 
derlying orbits and the input are somewhat different. Fig. [6] 
shows that the reconstructed solution at the minimum over- 
lies the input orbit almost exactly. We conclude that when 
the error bars are as small as in PD1, only one orbit is con- 
sistent with the data. 

Fig. [8] shows the results obtained with PD2, PD3 and 
PD4. For PD2 and PD3, which have small to moderate er- 
ror bars, the Metropolis algorithm reduces D' to the noise 
floor only for so in the range (44, 46) kpc. The scatter be- 
tween runs is small, cr^/ ~ 0.1. Fig. 1111 and Fig. 1131 show 
the reconstructed distances and tangential velocities asso- 
ciated with the best two solutions found at 44 and 46 kpc. 
With PD2, these reconstructions provide a distance esti- 
mate to the stream that is, at worst, 2 kpc in error, and a 
v± estimate that is at worst 5 kms -1 in error. With PD3, 
the reconstruction is, at worst, 3 kpc and 10 kms -1 in er- 
ror. For many sections of the orbits, the errors are less than 
stated. Thus the method can identify orbits consistent with 
the stream, and reject those that are inconsistent with it. 

For PD4, which has large error bars comparable to those 
for velocity data from the SDSS, scatter between repeat runs 
is much larger, <j D ' ~ 0.5. Only distances so < 40 kpc can be 
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Figure 8. The same as Fig. [5] but for input data sets PD2 (up- 
per), PD3 (middle) and PD4 (lower). In PD2 and PD3, the noise 
floor D' ~ —4 is reached only for starting distances in the range 
44 — 46 kpc. In PD4 we find an isolated good solution at 43 kpc, 
but in contrast to what happens with datasets PD2 and PD3, 
as we change distance the value of D' oscillates around ~ —3.25 
for most of the range, rather than varying smoothly. This be- 
haviour arises because the volume of parameter space that has 
to be searched is large on account of the largeness of the velocity 
errors. Such an extensive volume of parameter space cannot be 
exhaustively searched with the allocated computational resource. 
Only orbits with so < 40 kpc could be excluded with confidence. 



confidently excluded, although a high-quality solution near 
so ~ 43 kpc provides a reconstruction in error by at most 
2 kpc in distance and 5kms -1 in v±. Fig. 1111 and Fig. [13] 
shows us that ignoring this high-quality solution would give 
reconstructed quantities in error by up to 5 kpc in distance 
and 12kms _1 in v±. Increasing the errors associated with 
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Figure 9. The same as Fig.[8]but for input data sets PD5 (upper), 
PD6 (lower). The noise floor D ~ —3.5 is now approached over a 
wider range in so: 42 — 47 kpc in PD5 with some confidence, and 
38 — 48 kpc in PD6 with little confidence. As the errors increase, 
the search becomes a more arduous task, and patchy performance 
of this task is reflected in the rough bottoms to the graphs. 

the input clearly permits less satisfactory solutions to be 
returned, and also decreases the ability of the search proce- 
dure to find true orbits at particular distances, should they 
exist, by increasing the volume of parameter space it has to 
search. The former complaint is a physical statement about 
the limitations of the input data. The latter complaint is 
an algorithmic one, which may be remedied by providing 
the search with more computational cycles, or improving its 
efficiency. 

Fig. [9] shows the results obtained from input sets PD5 
and PD6 in which offsets were applied to the input velocities. 
The reconstructed distances and tangential velocities for in- 
teresting tracks are shown in Fig. [12] and Fig. [H For PD5, 
the scatter between runs is <jd> ~ 0.5, and the distance to 
the stream so can be said to lie in the range 42 — 47 kpc with 
some confidence. Fig. 1121 shows that the stream does indeed 
lie in this range, so the algorithm has successfully corrected 
the small offset. The error in reconstructed distance and v± 
would be 3 kpc and 8kms _1 at worst. 

PD6 demonstrates the limits of the method. On account 
of the large scatter in D' , a D > ~ 1, we cannot identify the 
correct distance from Fig. [9] We expect the distance range 
for permitted orbits to be wide with this input, because 
the velocity error bars are large. Fig. [9] illustrates another 
problem of using large error bars: the search becomes harder 
because the parameter space to be searched is much larger. 
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Figure 10. The same as Fig. but for input data set PD7. The 
results are very similar to those of PD2 and PD3 (Fig. [8}, with 
a noise floor D' ~ —4 and solutions acceptable only in the range 
44 - 46 kpc. 



Consequently the plots of the results have rough bottoms, 
where the algorithm has failed to reach consistent minima 
for searches at adjacent distances. This can be remedied 
by re-running the search with more deformations and more 
iterations, and may be addressed in future upgrades to the 
search procedure. 

PD7 demonstrates that the accuracy of the method is 
not necessarily significantly degraded when the number of 
velocity data points is substantially lower than the number 
of positional data points, as might be expected from a stream 
for which the only available radial velocities are those of 
giant stars. The results for PD7 are shown in Fig. [10] and 
are directly comparable to those of PD2 and PD3 (Fig. [5}. 
In particular, the range of allowed distances, 44 — 46 kpc, 
is the same and the reconstructed orbits show comparable 
accuracy (Figs. [12] and [T4|) . 



5 THE EFFECT OF CHANGING THE 
POTENTIAL 

Paper I demonstrated the ability of orbit reconstruction to 
diagnose the Galactic potential with astonishing precision 
when the track of an orbit on the sky is precisely known. 
Here we investigate the ability of orbit reconstruction to 
identify the correct potential when the track has to be in- 
ferred from realistic stream data. 

We use as our input the PD2 data set from Section [4] 
The top panel of Fig. [15] shows the results of asking the 
algorithm to find orbits in a Kepler potential with mass 
M = 4.18 x 10 11 Mq, which produces roughly the same pas- 
sage time along the stream as does Model II. The bottom 
panel shows the results obtained using a potential of the 
form $(r) = rf r , which gives a rotation curve of the form 
v c (r) = y/rfr, with f r — 6.86x 10 2 ( kms -1 ) 2 / kpc again cho- 
sen to produce the same passage time along the stream as 
does Model II. These two potentials represent, respectively, 
relatively extreme falling and rising rotation curve models. 
We do not offer them as realistic candidates for the Galactic 
potential, but intend to demonstrate that model potentials 
with approximately correct radial force, but incorrect shape, 
can be excluded using this method. 
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Figure 11. Heliocentric distance for selected reconstructed orbits 
from PD2, PD3 and PD4. The tracks selected are those with 
lowest D' at distances for which D' approaches the noise floor in 
Fig. ED 



The distance range considered in both cases spans ap- 
proximately ±15 percent of the true distance, which is the 
uncertainty one might expect in distances obtained by pho- 
tometry. Since all values of D' are several orders of mag- 
nitude larger than the values one obtains with the correct 
potential, it is clear that no orbits can be found at the dis- 
tances considered, and both potentials can be excluded. 

The upper panel of Fig. [16] shows the results of recon- 
structing an orbit from the PD2 data set in a potential that 
differs from the Model II potential used to define the tidal 
stream only in that the halo mass has been varied by the 
specified ratio. This high-quality data set yields a sharp min- 
imum in D' as a function of so (Fig. [8} and the upper panel 
of Fig. \W\ shows this minimum value of D' as a function 
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Figure 12. As Fig. [IT] except for data sets PD5, PD6 and PD7. 
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Figure 13. Tangential velocities for selected reconstructed orbits 
from PD2, PD3 and PD4. The tracks selected are the same as in 

Fig.HU 



of assumed halo mass. The minimum of D' lies at D' < — 4 
when the mass used lies within ~ 5 percent of the true value 
and D' > —3.85 otherwise. 

The lower panel of Fig. 1161 shows that the value of so at 
which D' attains its minimum for given halo mass decreases 
systematically as the halo mass increases. Insight into this 
behaviour can be obtained by considering the discretised 
equation of angular-momentum conservation 



A(sv±) 



d$ . 
-s- — At ■■ 
dr± 



d$ Av\\ 



(17) 



where A implies the change in a quantity between successive 
data points. By expanding the left side to first order in small 
quantities we can obtain an expression for Av±. Summing 
the changes in v± along the track we have 



, 1N , v ( d<F Av\ 

v ± (end) -v± (start) = - ^ ^^-pr 

An independent equation for v± is 



As 



v± 



Au 
'At 



: SF\\ 



Au 

Av\\ 



(18) 



(19) 



where we have used the radial equation of motion. Equa- 
tions (|18|) and (|19p yield independent estimates of v± (end) — 
v± (start). The right side of equation ([15)) yields an estimate 
that is independent of the scaling of $ but systematically 
decreases with increasing s, while the right side of equation 
([T9)) yields an estimate that is proportional to the scaling of 
<£, but is almost independent of s, since F\\ ocl/s. When the 
machine is asked to reconstruct the orbit with F\\ taken too 
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Figure 14. As Fig. [13] except for data sets PD5, PD6 and PD7. 
The tracks selected are the same as in Fig. ^] 



small, it can change the right side of equation (|18[) to match 
the new value of the right side of equation (|19[) by increas- 
ing s. In this way, the discrepancies between the left sides of 
equations (|18[) and (|19p . which contribute substantially to 
the diagnostic D' , can be largely eliminated by increasing s 
as $ is scaled down. 

In principle this variation in reconstructed distance with 
the scaling of the potential could be combined with photo- 
metric distances to constrain the potential. Unfortunately, 
Fig. [16] shows that in the particular geometry under consid- 
eration, even a 10 percent distance error would produce a 
~ 50 percent error in the estimate of the halo mass. Fur- 
ther work is required to discover what effect the geometry 
of a particular stream has on its sensitivity to the potential. 
Also of interest is whether simultaneously using multiple 
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Figure 15. As the top panel of Fig. [8] except the reconstruction 
(using PD2) takes place in a Keplerian potential (top panel), and 
(bottom panel) a potential with <E>(r) oc r. In both cases, the 
constants are set to generate approximately the same passage 
time along the stream as in Model II. Comparing with Fig. [8] 
shows the values of D' are very poor, demonstrating that no orbits 
can be found in these potentials, so the latter can be excluded. 

streams, or streams with multiple wraps (such as the Sagit- 
tarius Dwarf stream), can provide tight constraints on the 
potential. 



6 CONCLUSIONS 

Paper I demonstrated the tremendous diagnostic power that 
is available if one knows the track of an orbit on the sky 
and the associated line-of-sight velocities. Tidal streams are 
made up of objects that are on closely related orbits. In 
particular, they roughly delineate the underlying orbit, but 
they do not do so exactly. We have presented a technique 
for identifying an underlying orbit and thus predicting the 
dynamical quantities that have not been observed, such as 
the distances to the stream and the proper motions of its 
particles. 

The technique involves defining a space of tracks on the 
sky and sequences of line-of-sight velocities that are con- 
sistent with the observational data, given the observational 
errors and the extent to which streams deviate from orbits. 
The equations of Paper I are used to determine a candidate 
orbit for each track, and then the extent to which the can- 
didate satisfies the equations of motion is quantified - in 
Paper I only violations of energy conservation along a can- 
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Figure 16. The upper panel shows the minimum value of the 
diagnostic for attempted reconstruction of the PD2 input, versus 
the halo mass ratio of the Dehnen-Binney potential in which the 
reconstruction has been attempted. Mo here is the value of the 
halo mass in the Model II potential. The lower panel shows the 
characteristic distance, so, of the best solution, versus halo mass 
ratio. The error bars represent the approximate uncertainty in 
the recovered result. 



didate were quantified. This diagnostic quantity is then used 
to search for tracks that could be projections of orbits in the 
Galactic potential. In practice the search is conducted for 
several possible distances to a fiducial point on the stream. 
If constraints on this distance are available from photome- 
try, the computational effort of the search can be reduced 
by narrowing the range of distances for which searches need 
to be conducted. 

We have taken the Galaxy's gravitational potential and 
the solar velocity with respect to the Galactic centre to be 
known. Our test s revolve around an N -body model of the 
Orphan Stream (Belo kurov et al.ll2007l ). We show that for 
this stream, which is ~ 40° long and 0.3° wide at its ends, 
distances to and tangential velocities of points on the stream 
can be recovered to within ~ 2kpc and ~ 5kms -1 , respec- 
tively, if radial velocities accurate to ~ lkms -1 are mea- 
sured. As the errors in the measured radial velocities in- 
crease, the space of tracks that must be considered grows 
bigger and the search for acceptable orbits becomes more 
laborious. Moreover, the range of distances for which ac- 
ceptable orbits can be found broadens. However, even with 
errors in radial velocities as large as ±10kms _1 , the uncer- 
tainties in the recovered distances are no greater than ^10 
percent and the recovered tangential velocities are accurate 



to better than ~ 20 percent. Zero-point errors in the input 
velocities that are reflected in appropriately wide error bars 
broaden the range of acceptable orbits but do not skew the 
results. 

We have shown that the method maintains its accuracy 
even when very few radial velocity points are used to de- 
fine the input, as might be necessary when radial velocities 
can only be measured for giant stars. In our tests, compa- 
rable results were obtained from pseudodata based on only 
three velocity measurements and from pseudodata based on 
fifteen velocity measurements. Indeed, the results obtained 
with three accurate velocity measurements were significantly 
superior to those obtained with fifteen lower-quality mea- 
surements. Naturally, exactly how many points are required 
to provide well-determined input will depend on the shape 
of the radial velocity curve along the stream in question. 

We expect the accuracy of reconstructions to depend 
on the geometry of the problem in hand. In particular, we 
expect streams at apocentre, where families of orbits are 
compressed both on the sky and in radial velocity, to yield 
poorer results than streams away from apocentre. Unfortu- 
nately, streams are most likely to be discovered at apocentre 
because both orbital compression and low proper motions 
around apocentre lead to a high density of stars at apoc- 
entre. We expect streams that are relatively narrow to pro- 
duce more accurate results, because the permitted deviation 
of the orbit from the stream is then low. We also expect to 
have more difficulty reconstructing orbits from streams that 
contain a visible progenitor, since the potential of the pro- 
genitor will cause orbits in the progenitor's vicinity to differ 
materially from orbits in the Galaxy's underlying potential. 

Paper I suggested that it should be possible, if suffi- 
ciently accurate input is provided, to constrain the Galac- 
tic potential, since the wrong potential will not admit an 
acceptable orbit. We have tested this possibility for input 
with realistic errors. We find that two potentials of signifi- 
cantly different shape, the Kepler potential and $(r) oc r, 
are clearly excluded. We have also tested for changes in scal- 
ing of an otherwise correctly-shaped potential, by varying 
the mass of the assumed dark halo around the value used 
to make the pseudo-data. In this case, we find the correct 
potential is identified, with the diagnostic quantity gener- 
ally worsening as the halo mass moves away from its correct 
value by more than ~ 5 percent. We further find a consis- 
tent relationship between the reported stream distance and 
the halo mass with which the reconstruction takes place. 
Although the reported distance is only weakly dependent 
upon halo mass, this does open the possibility of using al- 
ternative distance measurements, such as photometric dis- 
tances, in conjunction with these techniques to constrain 
the Galactic potential. Further work is necessary to deter- 
mine a full scheme to recover parameters of the potential 
from stream data. Also in question is the extent to which 
simultaneous reconstruction of multiple streams, and recon- 
struction of streams with multiple wraps around the Galaxy, 
might provide stronger constraints on the Galactic potential 
than the short section of a single wrap that we have consid- 
ered. It may also prove possible to refine the other main 
assumption of our scheme, the location and velocity of the 
Sun. 

It is instructive to compare our method of finding orbits 
of progenitors with the traditional N-body method. First our 
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method explores each orbit at a tiny fraction of the compu- 
tational expense of N-body modelling, so it is feasible to 
automate the search of orbit space. Moreover, the search 
lends itself to parallelisation. Whereas only a successful at- 
tempt to model a stream with N bodies yields an interesting 
conclusion, our method can show that no orbit is consistent 
with a given range of distances. 

There are several directions in which this work could be 
profitably extended. 

1. There is scope for a powerful synergy between our 
method and N-body modelling: our method is first used 
to identify a likely orbit and this orbit then provides ini- 
tial conditions for an N-body simulation, which reveals the 
offset between the progenitor's orbit and the stream. This 
knowledge would enable the bow-tie region to be made nar- 
rower. Finally our method is used again to determine the 
orbit with still higher precision. 

2. Currently VLBI observations of masers yield trigono- 
metric parallaxes for of order two dozen sources located at 
several kpc from t he Sun that are accurate to several percent 
(jReid et al.ll20Q9h . and Gaia will yield results of similar pre- 
cision for several million stars. The method discussed here 
promises distances of slightly higher precision to sources that 
lie 50 to 100 kpc from the Sun. T hese "geometrodynamicaF 
distances in the terminology of iJin &; Lvnden-Belll (2007) 
may over time play a significant role in astrophysics, just 
as trigonometric distances did before them, by checking and 
calibrating photometric distances. However, before this ex- 
citing prospect can be realised we must overcome the prob- 
lem that the method merely links distances, velocities and 
the still uncertain gravitational potential of the Galaxy. How 
can we most effectively exploit this link between measure- 
ments of radial velocities, proper motions and photometry to 
obtain tight constraints on both the distances and the poten- 
tial? This is an extremely important question given current 
interest in mapping the Galaxy's dark matter through the 
gravitational force field that it generates. Undoubtedly, pin- 
ning down the potential will be greatly facilitated by mod- 
elling several streams simultaneously. 

3. Finer discrimination between candidate orbits would be 
possible if one could lower the noise floor on the diagnostic 
function D' by upgrading the numerical methods used to 
obtain solutions of equations (pQ). The scheme for searching 
the space of possible tracks could also be made faster and 
more reliable. 

The list of streams to which this technique could be 
applied is already quite long: obvious examples include the 
tidal tails of the globular clusters Palomar 5 and NGC 5466, 
those of the Large Magellanic Cloud, the Orphan Stream, 
and the tidal stream of the Sagittarius dwarf galaxy. How- 
ever, before we can exploit the methods presented here, ra- 
dial velocities are needed along these streams. The greatest 
precision is promised by the narrowest streams, and these 
are well defined only in main-sequence stars. Therefore the 
observational challenge is to obtain accurate radial veloci- 
ties for tens of faint stars. This will require 8 m telescope 
time for what may seem unglamorous work. This paper sug- 
gests, however, that the scientific rewards of such observa- 
tions would be far-reaching. 

Recently it has been realised that orbits reconstruction 
is possible using proper motions along the stream rather 



than radial velocities (|Evre & Binnevll2009h . The principles 
developed in the present paper will undoubtedly transfer 
to orbit reconstructions from proper motions. It remains to 
be discovered whether they will enable useful distances to 
be determined from proper motions of currently attainable 
precision. We hope to report on this issue shortly. 
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